Constraining neutron star tidal Love numbers with gravitational wave detectors 
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Ground-based gravitational wave detectors may be able to constrain the nuclear equation of state 
using the early, low frequency portion of the signal of detected neutron star - neutron star inspirals. 
In this early adiabatic regime, the influence of a neutron star's internal structure on the phase of the 
waveform depends only on a single parameter A of the star related to its tidal Love number, namely 
the ratio of the induced quadrupole moment to the perturbing tidal gravitational field. We analyze 
the information obtainable from gravitational wave frequencies smaller than a cutoff frequency of 
400 Hz, where corrections to the internal-structure signal are less than 10%. For an inspiral of two 
non-spinning 1.4M© neutron stars at a distance of 50 Megaparsecs, LIGO II detectors will be able to 
constrain A to A ^ 2.0 x 10 37 gcm 2 s 2 with 90% confidence. Fully relativistic stellar models show that 
the corresponding constraint on radius R for 1.4M© neutron stars would be R ^ 13.6 km (15.3 km) 
for a n = 0.5 (n = 1.0) polytrope with equation of state p oc p 1+1//n . 



Background and motivation: Coalescing binary neutron 
stars are one of the most important sources for gravi- 
tational wave (GW) detectors fl. LIGO I observations 
have established upper limits on the event rate , and at 
design sensitivity LIGO II is expected to detect inspirals 
at a rate of ~ 2/day 

One of the key scientific goals of detecting neutron 
star (NS) binaries is to obtain information about the 
nuclear equation of state (EoS), which is at present 
fairly unconstrained in the relevant density range p ~ 
2 — 8 x 10 14 gcm -3 The conventional view has been 
that for most of the inspiral, finite-size effects have a 
negligible influence on the GW signal, and that only dur- 
ing the last several orbits and merger at GW frequencies 
/ ^ 500 Hz can the effect of the internal structure be 
seen. 

There have been many investigations of how well the 
EoS can be constrained using these last several orbits and 
merger, including constraints from the GW energy spec- 
trum 5], and, for black hole/NS inspirals, from the NS 
tidal disruption signal [6| . Several numerical simulations 
have studied the dependence of the GW spectrum on the 
radius 7]. However, there are a number of difficulties as- 
sociated with trying to extract equation of state informa- 
tion from this late time regime: (i) The highly complex 
behavior requires solving the full nonlinear equations of 
general relativity together with relativistic hydrodynam- 
ics, (ii) The signal depends on unknown quantities such 
as the spins of the stars, (iii) The signals from the hydro- 
dynamic merger (at frequencies > 1000 Hz) are outside 
of LIGO's most sensitive band. 

The purpose of this paper is to demonstrate the po- 
tential feasibility of instead obtaining EoS information 
from the early, low frequency part of the signal. Here, 
the influence of tidal effects is a small correction to the 
waveform's phase, but it is very clean and depends only 
on one parameter of the NS - its Love number [8| . 

Tidal interactions in compact binaries: The influence of 
tidal interactions on the waveform's phase has been stud- 



ied previously using various approaches @, 0, E3, eh, m 

Il3j |. We extend those studies by (i) computing the ef- 
fect of the tidal interactions for fully relativistic neutron 
stars, i.e. to all orders in the strength of internal gravity 
in each star, (ii) computing the phase shift analytically 
without the assumption that the mode frequency is much 
larger that the orbital frequency, and (iii) performing a 
computation of how accurately the Love number can be 
measured. 

The basic physical effect is the following: the I = 2 
fundamental /-modes of the star can be treated as forced, 
damped harmonic oscillators driven by the tidal field of 
the companion at frequencies below their resonant fre- 
quencies. Assuming circular orbits they obey equations 
of motion of the form 20J1 



q + 7? + Lolq = A(t) cos[m$(t)] 



(1) 



where q(t) is the mode amplitude, loq the mode frequency, 
7 a damping constant, m is the mode azimuthal quan- 
tum number, <fr(t) is the orbital phase of the binary, 
and A(t) is a slowly varying amplitude. The orbital fre- 
quency Lo(t) — $ and A(t) evolve on the radiation reac- 
tion timescale which is much longer than 1/u)q. In this 
limit the oscillator evolves adiabatically, always tracking 
the minimum of its time-dependent potential. The en- 
ergy absorbed by the oscillator up to time t is 
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The second term here describes a cumulative, dissipative 
effect which dominates over the first term for tidal in- 
teractions of main sequence stars. For NS-NS binaries, 
however, this term is unimportant due to the small vis- 
cosity 11[, and the first, instantaneous term dominates. 

The instantaneous effect is somewhat larger than of- 
ten estimated for several reasons: (i) The GWs from the 
time varying stellar quadrupole are phase coherent with 
the orbital GWs, and thus there is a contribution to the 
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energy flux that is linear in the mode amplitude. This af- 
fects the rate of inspiral and gives a correction of the same 
order as the energy absorbed by the mode ■ (ii) Some 
papers [t| [ll], [l2| compute the orbital phase error as a 
function of orbital radius r. This is insufficient as one has 
to express it in the end as a function of the observable fre- 
quency, and there is a correction to the radius-frequency 
relation which comes in at the same order, (iii) The ef- 
fect scales as the fifth power of neutron star radius R, 
and most previous estimates took R = 10 km. Larger NS 
models with e.g. R — 16 km give an effect that is larger 
by a factor of ~ 10. 

Tidal Love number: Consider a static, spherically sym- 
metric star of mass m placed in a time-independent ex- 
ternal quadrupolar tidal field Eij. The star will develop 
in response a quadrupole moment Qij. In the star's lo- 
cal asymptotic rest frame (asymptotically mass centered 
Cartesian coordinates) at large r the metric coefficient 
gu is given by (in units with G — c = 1) [141 ] : 
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(3) 



where n % — x l /r; this expansion defines the traceless ten- 
sors Eij and Qij. To linear order, the induced quadrupole 
will be of the form 
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Here A is a constant which we will call the tidal Love 
number (although that name is usually reserved for the 
dimensionless quantity k% = |GAi?~ 5 ). The relation 
(HI between and £^ defines the Love number A for 
both Newtonian and relativistic stars. For a Newto- 
nian star, (1 — g u ) /2 is the Newtonian potential, and 
Qij is related to the density perturbation dp by Qij = 
J d 3 x5p {xiXj — r 2 5ij/3). 

We have calculated the Love numbers for a variety of 
fully relativistic NS models with a polytropic pressure- 
density relation P = Kp 1+1 / n . Most realistic EoS's re- 
semble a poly trope with effective index in the range n ~ 
0.5 — 1.0 [15| . The equilibrium stellar model is obtained 
by numerical integration of the Tolman-Oppcnhcimer- 
Volkhov equations. We calculate the linear I — 2 static 
perturbations to the Schwarzschild spacetimc following 
the method of The perturbed Einstein equations 

5G^ — % , k§T ix u can be combined into a second order dif- 
ferential equation for the perturbation to gu- Matching 
the exterior solution and its derivative to the asymptotic 
expansion ([3]) gives the Love number. For m/R ~ 10~ 5 
our results agree well with the Newtonian results of Refs. 
@; [I3|- Figure 1 shows the range of Love numbers for 
m/R = 0.2256, corresponding to the surface redshift 
z = 0.35 that has been measured for EXO0748-676 [Hj]. 
Details of this computation will be published elsewhere. 

Effect on gravitational wave signal: Consider a binary 
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FIG. 1: [Top] The solid lines bracket the range of Love num- 
bers A for fully relativistic polytropic neutron star models 
of mass m with surface redshift z = 0.35, assuming a range 
of 0.3 < n < 1.2 for the adiabatic index n. The top scale 
gives the radius R for these relativistic models. The dashed 
lines are corresponding Newtonian values for stars of radius 
R. [Bottom] Upper bound (horizontal line) on the weighted 
average A of the two Love numbers obtainable with LIGO II 
for a binary inspiral signal at distance of 50 Mpc, for two 
non-spinning, 1.4Mq neutron stars, using only data in the 
frequency band / < 400 Hz. The curved lines are the actual 
values of A for relativistic polytropes with n = 0.5 (dashed 
line) and n = 1.0 (solid line). 



with masses mi, mi and Love numbers Ai, A2. For sim- 
plicity, we compute only the excitation of star 1; the 
signals from the two stars simply add in the phase. Let 
w„, Ai,n and Qfj be the frequency, the contribution to 
Ai and the contribution to Qij of modes of the star with 
I = 2 and with n radial nodes, so that Ai = S n Ai jn 
and Qij — S„Q™ . Writing the relative displacement as 
x = (r cos $, r sin $,0), the action for the system is 
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(5) 



Here M and /x are the total and reduced masses, and 
Eij = —m-ididj (l/r) is the tidal field. This action is 
valid to leading order in the orbital potential but to all 
orders in the internal potentials of the NSs, except that it 
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neglects GW dissipation, because Qij and _£h are denned 
in the star's local asymptotic rest frame [19( . 

Using the action ((SJ), adding the leading order, Burke- 
Thorne GW dissipation terms, and defining the total 
quadrupole Qjj — + iiXiXj — [ir 2 5ij/3 with = 
£„<2™ , gives the equations of motion 
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(6a) 



1 2 d 5 Qj- 
m 2 X 1 . n uj 2 l d i d j - - ~\i,„uJn 5 i3 .(6b) 



By repeatedly differentiating and eliminating second 
order time derivative terms using the conservative parts 
of Eqs. ©, we can express d 5 Qjj/dt 5 in terms of x % , if, 
Q^j and and obtain a second order set of equations; 
this casts Eqs. ([5]) into a numerically integrable form. 

When GW damping is neglected, there exist equi- 
librium solutions with r = const, $ = $o + wt for 
which Qjj is static in the rotating frame. Working to 
leading order in Ai, n , we have = Q! + Qcos(2$), 
gT = Q' - Qcos(2<I>), Qj 2 = Qsin(20>), Qj 3 = -2Q', 
where 



Q = \^ 



3™ 2 Ai,„ , _ 1 2 
2(l-4a; 2 )r 3 ' J ~ 6^ 
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(7) 

and x n — u>/uj n . Substituting these solutions back 
into the action ([5]), and into the quadrupole formula 

E = —%{QijQij) for the GW damping, provides an effec- 
tive description of the orbital dynamics for quasicircular 
inspirals in the adiabatic limit. We obtain for the orbital 
radius, energy and energy time derivative 
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where \n — m 2 \i_ n uj w /' i m 1 1 M 5 / 3 , gx(x) = 1 + 
3/(1 - 4a; 2 ), g 2 (x) = 1 + (3 - 4x 2 )(l - 4a; 2 )" 2 , and 
g 3 (x) = (M/m.2 + 2 - 2a; 2 ) /(l - 4a; 2 ). Using the for- 
mula d 2 ^/duj 2 = 2 {dE/duj) /E for the phase *(/) of the 
Fourier transform of the GW signal at GW frequency 
/ = lo/tt |25j now gives for the tidal phase correction 

W) = -l^T, X ^l V dv'v'(v^v^ gi (x' n ) 1 

54(35) 
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(9) 
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FIG. 2: [Fop] Analytic approximation (|10p to the tidal pertur- 
bation to the gravitational wave phase for two identical IAMq 
neutron stars of radius R = 15 km, modeled as n = 1.0 poly- 
tropes, as a function of gravitational wave frequency /. [Bot- 
tom] A comparison of different approximations to the tidal 
phase perturbation: the numerical solution (lower dashed, 
green curve) to the system ©, and the adiabatic analytic 
approximation @ (upper dashed, blue), both in the limit 
(|11[) and divided by the leading order approximation (|10f) . 



Here v = (nMf) 1 / 3 , Vi is an arbitrary constant related 
to the initial time and phase of the waveform, and x' n = 
(v 1 ) 3 1 (Mu) n ). In the limit uj <C uj n assumed in most 
previous analyses @, @, [II], [H} , we get 
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(10) 



which depends on internal structure only through Ai and 
\ 2 . Here we have added the contribution from star 2. 
The phase (TT0|) is formally of post-5-Newtonian (P5N) 
order, but it is larger than the point-particle P5N terms 
(which are currently unknown) by ~ (R/M) 5 ~ 10 5 . 

Accuracy of Model: We will analyze the information con- 
tained in the portion of the signal before / = 400 Hz. 
This frequency was chosen to be at least 20% smaller 
than the frequency of the innermost stable circular orbit 
[ill ] for a conservatively large polytropic NS model with 
n = 1.0, AI = 1.4M Q , and R = 19km. We now argue 
that in this frequency band, the simple model (jTUJ) of the 
phase correction is sufficiently accurate for our purposes. 

We consider six types of corrections to (flQ|) . For 
each correction, we estimate its numerical value at the 
frequency / = 400 Hz for a binary of two identical 
m = 1.4M Q , R = 15, n = 1.0 stars: (i) Corrections 
due to modes with I > 3 which are excited by higher or- 
der tidal tensors £yfc, .... The I — 3 correction to E(u>), 
computed using the above methods in the low frequency 
limit, is smaller than the I = 2 contribution by a factor 
of 65k 3 R 2 /(45k 2 r 2 ), where k 2 , k% are apsidal constants. 
For Newtonian polytropes we have k 2 = 0.26, k% = 0.106 
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Q, and the ratio is 0.58(i?/r) 2 = 0.04(i?/15 km) 2 , (ii) 
To assess the accuracy of the w < w n limit underlying 
(fTTTj) wc simplify the model §5$ by taking 



The rms statistical measurement 



for all n, 



(11) 



so that Qfj/Xx !n is independent of n. This simplification 



does not affect ([TO]) and increases the size of the finite fre- 
quency corrections in (|9|) since uj n > luq 1 . This will yield 
an upper bound on the size of the corrections. (Also the 
n > 1 modes contribute typically less than 1 — 2% of the 
Love number § .) Figure 2 shows the phase correction 6^ 
computed numerically from Eqs. (O, and the approxima- 
tions flH]) and (TTU)) in the limit (fTTj) . We see that the adia- 
batic approximation @ is extremely accurate, to better 
than 1%, and so the dominant error is the difference be- 
tween © and (jTOjl . The fractional correction to (jTOj) 
is ~ 0.7x 2 ~ 0.2(///o) 2 , where /„ = u /(2w), neglect- 
ing unobservable terms of the form a + 0f. This ratio is 
< 0.03 for / < 400 Hz and for f a > 1000 Hz as is the case 
for /-mode frequencies for most NS models ■ (hi) We 
have linearized in Ai; the corresponding fractional cor- 
rections scale as {R/rf ~ 10" 3 (i?/15 km) 5 at 400 Hz. 

(iv) The leading nonlinear hydrodynamic corrections can 
be computed by adding a term — aQ^Q^Q^ to the La- 
grangian ([5]), where a is a constant. This corrects the 
phase shift JTOJ) by a factor l-285aA 2 nwV968 ~ 0.9995, 
where we have used the models of Ref. [21| to estimate a. 

(v) Fractional corrections to the tidal signal due to spin 
scale as ~ / S 2 pin //L. where / spin is the spin frequency 
and / m ax the maximum allowed spin frequency. These 
can be neglected as / max > 1000 Hz for most models 
and / S pin is expected to be much smaller than this, (vi) 
Post-l-Newtonian corrections to the tidal signal (fT0|) will 
be of order ~ M/r ~ 0.05. However these corrections 
will depend only on Ai when oj «C u) n , and can easily be 
computed and included in the data analysis method we 
suggest here. 

Thus, systematic errors in the measured value of A due 
to errors in the model should be < 10%, which is small 
compared to the current uncertainty in A (see Fig. 1). 

Measuring the Love Number: The binary's parameters 
are extracted from the noisy GW signal by integrating 
the waveform h(t) against theoretical inspiral templates 
h(t,9' 1 ), where 8 l are the parameters of the binary. The 
best-fit parameters 9 l are those that maximize the over- 
lap integral. The probability distribution for the sig- 
nal parameters for strong signals and Gaussian detec- 
tor noise is p (0 l ) = A/exp (-1/2 r^-AcMc^) where 
A8 l = 6 l - 9\ Y i} = (dh/dd* , dh/d6i ) is the Fisher in- 
formation matrix, and the inner product is defined by 



1 Buoyancy forces and associated g-modes for which u) n < uio have 
a negligible influence on the waveform's phase|2d(|. 



Eq. (2.4) of Ref. [2j 

error in 9 l is then *J (r _1 ) M . 

Using the stationary phase approximation and neglect- 
ing corrections to the amplitude, the Fourier transform 
of the waveform for spinning point masses is given by 
h(f) = „4f~ 7/6 exp (£*)■ Here the phase * is 
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where v = (7rM/) 1//3 , and a are spin parameters, and 
7 is Euler's constant [22j. The tidal term (TIT)]) adds lin- 
early to this, yielding a phase model with 7 parameters 
(t c ,(j) c ,M,fi,0,(T,X), where A = [(llm 2 + M)X\/mi + 
(llmi + M)A 2 /m 2 ]/26 is a weighted average of \i and 
A 2 . We incorporate the maximum spin constraint for 
the NSs by assuming a Gaussian prior for and a as in 
Ref. [13. 

Figure 1 [bottom panel] shows the 90% confidence up- 
per limit A < 20.1 x 10 36 gcm 2 s 2 we obtain for LIGO 
II (horizontal line) for two nonspinning 1.4M Q NSs at a 
distance of 50 Mpc (signal-to-noise of 95 in the frequency 
range 20 — 400Hz) with cutoff frequency f c = 400 Hz, as 
well as the corresponding values of A for relativistic poly- 
tropes with rt = 0.5 (dashed curve) and n = 1.0 (solid 
line). The corresponding constraint on radius assuming 
identical 1.4M Q stars would be R ^ 13.6 km (15.3 km) 
for n = 0.5 (n = 1.0) polytropes. Current NS models 
span the range 10 km < R < 15 km. 

Our phasing model (JT2J) is the most accurate avail- 
able model, containing terms up to post-3. 5-Newtonian 
(P3.5N) order. We have experimented with using lower 
order phase models (P2N, P2.5N, P3N), and we find that 
the resulting upper bound on A varies by factors of order 
~ 2. Thus there is some associated systematic uncer- 
tainty in our result. To be conservative, we have adopted 
the most pessimistic (largest) upper bound on A, which 
is that obtained from the P3.5N waveform. 

In conclusion, even if the internal structure signal is 
too small to be seen, the analysis method suggested here 
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could start to give interesting constraints on NS internal 
structure for nearby events. 
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